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The dynamics of a discrete soliton in an array of Bose-Einstein condensates under the action of 
a periodically time-modulated atomic scattering length ('Feshbach-resonance management, FRM") 
is investigated. The cases of both slow and rapid modulation, in comparison with the tunneling 
frequency, are considered. We employ a discrete variational approach for the analysis of the system. 
The existence of nonlinear resonances and chaos is predicted at special values of the driving fre- 
quency. Soliton splitting is observed in numerical simulations. In the case of the rapid modulation, 
we derive an averaged equation, which is a generalized discrete nonlinear Schrodinger equation, in- 
cluding higher-order effective nonlinearities and intersite nonlinear interactions. Thus the predicted 
discrete FRM solitons are a direct matter-wave analog of recently investigated discrete diffraction- 
managed optical solitons. 
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I. INTRODUCTION 

Discrete solitons in nonlinear lattices with periodically varying parameters have recently attracted much attention. 
System belonging to this type include arrays of optical waveguides subject to periodic diffraction management 1 ' 2 
and waveguide arrays with a periodic modulation of the tunnel-coupling constant 3 . The corresponding model is 
typically based on the discrete nonlinear Schrodinger (DNLS) equation, with the coefficient in front of the second- 
finite-difference term varying along the propagation distance (formally, it looks like periodic time modulation of the 
coefficient). It was shown that, in the case of rapid and strong variations of the coupling constant, a stable breathing 
discrete soliton can exists 2 (the so-called diffraction-managed soliton) . On the other hand, application of a relatively 
slow weak or moderate modulation at a resonant frequency results in a splitting of the discrete soliton 3 . 

In periodically modulated DNLS systems of another type, the coefficient of the on-site cubic nonlincarity is subject 
to the modulation. In terms of nonlinear optics, these may be arrays of waveguides which have a layered structure, 
with the strength 4 , or even sign 5 , of the nonlinearity alternating between layers. An alternative, and actually more 
straightforward, physical realization of this type of the lattice is offered by an array of droplets of a Bose-Einstein 
condensate (BEC) trapped in a deep optical lattice 6 ' 7 , with the BEC scattering length oscillating in time. The lat- 
ter type of the time-modulation may be provided by ac magnetic field tuned to the Feshbach resonance, as it was 
predicted theoretically 8 and demonstrated experimentally 9 . By analogy with the well-known techniques of the dis- 
persion management 10 and the above-mentioned diffraction management 1 in nonlinear optics, this time-modulation 
technique, applied to BEC, may be called Feshbach-resonance management (FRM). Very recently, it has been demon- 
strated that FRM provides for an effective mechanism of stabilization of two-dimensional BECs, even in the absence 
of the dc-magnetic-field trap 11 . One-dimensional solitons subject to the action of FRM were also recently studied, 
which reveals stable breathers oscillating between the Gaussian and Thomas-Fermi shapes, and stable breathers of 
other types 12 . 

The aim of the present work is to consider the dynamics of solitons in the one-dimensional DNLS model with the 
nonlinearity subject to periodic modulation. We will treat the cases of both relatively slow and rapid modulations. In 
the former case, we will apply an analytical variational approximation (VA), which was developed for one-dimensional 
lattice models in Ref. 13 , and direct simulations, to study resonances and splitting in the discrete-soliton dynamics 
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(a recent review of the VA technique can be found in Ref. 14 ). In the latter case, using the multiscale method 15 , we 
will derive an averaged equation, which has the form of a generalized DNLS equation with new nonlinear on-site and 
inter-site terms. Using this equation and VA, we will analyze the structure of average discrete-soliton solutions. 



II. THE MODEL 



We formulate the model in terms of a BEC trapped in a deep optical lattice, which is created by the interference 
of two counterpropagating optical beams. The dynamics of a BEC is governed by the Gross-Pitaevskii equation 16 : 

dib h 2 

M-£ = -^i> + v(TW + g (tM 2 i> , (1) 

where V(r) = Vq(x, y) sin 2 (fcz) is the optical potential, and g{t) — ^irh 2 a s {t) jm. Here a s is the time-varying atomic 
scattering length, and m is the atomic mass. As it was mentioned above, the time dependence of a s can be induced 
by ac magnetic field (or laser radiation with a time-modulated intensity) applied to the condensate. Due to the 
periodicity of V(r), for a weakly coupled array of BECs, one can present a solution as 

^ = 2«n(t)*(r-r„) , (2) 

n 

where the function <i>(r — r„) is assumed to be strongly localized around n-th site. Substituting Eq. (2) into Eq. (1), 
integrating over the transverse coordinates, and taking into account the exchange integrals only for neighboring sites, 
one arrives at a DNLS equation with a variable coefficient in front of the nonlinear term 6 ' 7 : 

1 2 
iu n + - (u„+i + tt n _i - 2u n ) + a(t) \u n \ u n = . (3) 

Here the overdot stands for the time derivative, time is made dimensionless by means of the rescaling t — > th/K, 
where K is the tunnel- coupling parameter between adjacent wells in the optical lattice 6,7 , and 

a(t) = a + di sin (ujt) (4) 

is a coefficient proportional to minus the atomic scattering length in the BEC. Equation (3) describes the dynamics 
restricted within the lowest Bloch zone. Account of inter-zone transitions requires an extension of the DNLS model 17 . 

Though the BEC system described above is the most relevant physical realization of Eq. (3), the same model also 
applies to an array of periodically modulated optical waveguides, with t being the propagation distance, rather than 
time. Without loss of generality, one can set a = 1 and a = — 1 in Eq. (4) for the cases of the negative and positive 
scattering lengths (repulsion and attraction between atoms), respectively. The wave function u n (t) is normalized so 
that the dynamical invariant of Eq. (3), 

oo 

W= Yl M 2 ' ( 5 ) 

n— — oo 

is the total number of particles. The characteristic length of the system 2ir/k ~ 1 /jm, Vq rs h 2 k 2 /m, and the 
atomic population in each well is ~ 10 3 atoms. The characteristic frequency for the tunneling between wells is 
CIl = K/h <~ 10 Hz, and the separation between the energy levels in a single well is ft > 100 Hz. Therefore, it makes 
sense to consider the variation of the driving frequency w in the interval Ql < uK/h < 0. 

First, we consider stationary pulse-shaped solutions of the unperturbed DNLS equation with a\ = in Eq. (4) in 
the form (see, e.g., Ref. 18 ' 19 ) 

u„(t) = Q n exp(mn - i\t) , (6) 

where k is a wave number and \ is a frequency. As it follows from the dispersion relation for the linearized equation (3), 
a localized solution with the maximum of Q n centered at some fixed point exists only for particular values of k at 
which the group velocity vanishes, so we take k = for ao = 1, or, equivalently, k = it for ao = — 1. 

The fundamental soliton for ao = 1 was studied in detail as a numerical solution to the nonlinear eigenvalue problem 
with zero boundary conditions 18-20 (see also a review in Ref. 21 ). In the case of ao = — 1, solitons are staggered 20 , with 
the 7r phase difference between adjacent sites; thus, on the contrary to the continuum nonlinear Schrodingcr (NLS) 
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equation, the DNLS model supports stable bright solitons for either sign (repulsion or attraction) of the nonlinear 
interaction. 

For convenience, here we briefly recapitulate basic properties of DNLS solitons. Parameters of the discrete soliton 
(6), found numerically from the nonlinear eigenvalue problem 18 , are shown by points in Fig. 1. The left axis in 
Fig. 1(a) pertains to the inverse width a, which was found by matching the soliton's tail to the asymptotic expression 
\u n \ = Acxp(— a|n|), where A is the soliton's amplitude [cf. Eq. (9) below]. The right axis in Fig. 1(b) corresponds 
to the pulse's area, which we define as S — J2 n \ Un \- The dependencies for ao = — 1 have the same form, with the 
only difference that \ IS shifted so that \ ~~ * 2 — \- Since W, \ an d ct are monotonic functions of A, the stationary 
solution (6) is defined by fixing of any one of these parameters. 

Similarly to the case of the continuum NLS equation, the addition of chirp to the soliton (6) (chirp imprinting) splits 
it into two separating soliton-like pulses, if the chirp b exceeds a critical (threshold) value 6 th (detailed consideration 
of a similar problem in the continuum NLS equation was given in Ref. 22 ). We introduce the chirp by taking an initial 
condition as 

u n (0) = Q n exp(ib\n\) (7) 

(the value of b is restricted to the interval [— tt,tt]). The dependence of 6 t h on the amplitude A of the unperturbed 
DNLS soliton is presented in Fig. 2. 

In fact, the curves shown in Fig. 2 diverge at sufficiently large A. The meaning of this is that, if the soliton's 
amplitude A exceeds the value 1.66, the initial pulse with any amount of chirp gives rise to a soliton centered at 
n = 0, while other parts of the initial pulse split off from it and move in opposite directions. 

It is possible to understand the chirp-induced splitting of the pulse into two in the following way. The original 
chirped pulse (7) may be regarded as a superposition of two pulses which carry the phase gradient of opposite signs 
(cf. a similar model developed in the framework of the continuum NLS equation in Ref. 22 ). As is known, the velocity 
of an isolated soliton is generated by its phase gradient. Since the two constituents of the overall chirped pulse are 
originally close to each other, the attraction between them is strong enough to keep them together. However, the 
increase of b leads to increase of the opposite phase-gradient thrusts applied to the constituents, and finally to splitting 
between them. 



III. THE VARIATIONAL APPROXIMATION AND DIRECT SIMULATIONS IN THE CASE OF SLOW 

MODULATIONS 

A. The general formalism of the variational approximation 

The DNLS equation (3) is derived from the Lagrangian, 

oo 
n— — oo 

^\u n+1 - u n \ 2 + ]^a(t)\u n \ A . (8) 

Following Ref. 13 , we base the variational approximation (VA) for the soliton governed by Eq. (3) on the following 
ansatz: 

u n (t) = Acxp(i<j) + ib\n\ — a\n\) , (9) 

where A, <f), b, and a are real functions of time. Substituting the ansatz (9) into Eq. (8), one can easily calculate the 
corresponding effective Lagrangian in a form 

L 1 db cos b 1 r . , sinha , . , . . 

777 = -7-77-T-T + i — + lWa(t) 3— cosh (2a , (10) 

W sinh(2a) dt cosha 4 cosh 3 a 

where 

IT = ^ 2 cotha, (11) 

is a dynamical invariant, which coincides with the total number of particles, obtained by substitution of the ansatz 
(9) into Eq. (5). We mention that a term in the full Lagrangian from which it follows that dW/dt — contains the 
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phase derivative <j> [which gives the frequency — \ in the stationary state, see Eq. (6)]. That term was dropped in 

the expression (10), as it does not contribute to other variational equations. Finally, the variational equations for the 
soliton's chirp and inverse width are 

db . , , sinh 3 a 1 r , ., l2 \ 2 cosh (2a) — 1 , . 

— = 2 (cos 6 — — — - -Wait) (tanh 2 a) — , (12) 

dt K ' cosh (2a) 2 w 1 ' cosh (2a) ' v ' 

dcy. 

— = - (sin b) (sinh a) tanh (2a) . (13) 



B. Revisiting the stationary model 

First, we dwell on the unperturbed case, with a(t) = const = a [cf. Eq. (4)]. In the case, all the points with a = 
and b = const are stationary solutions, i.e., fixed points (FPs). However, they do not correspond to localized waves, 
therefore they are formal solutions. Further, it is easy to see that Eqs. (12) and (13) give rise to nontrivial FPs with 
bpp = for ao = 1, and bpp = ir for do = — 1, the corresponding value app being defined by the equation 13 

sinha FP = (l + 3 tanh 2 a FP ) . (14) 

The parameters of the FP, which corresponds to stationary discrete soliton (6), are shown by solid lines in Fig. 1. 
As is seen, the results of the VA are in good agreement with the exact numerical solution of Eq. (3). Deviation in 
S (area of the pulse) indicates that the VA is not applicable in the limit of small A. This is clear because this limit 
corresponds to the continuum system, whose stationary soliton solution differs from the ansatz (9). 
Linearization of Eqs. (12) and (13) around the FP yields a squared frequency of small oscillations, 

2 sinh 3 (app ) cosh 2 (app ) 
W ° ~~ cosh 3 (2a FP ) 

x {4sinh(ap P )[cosh(2aFp) + 2] 

^ [5 cosh 2 (2a FP ) - 2 cosh(2a FP ) - l] 1 . (15) 

cosh (app) J 

Using Eq. (14), one can show that w 2 given by Eq. (15) is always positive, i.e., VA does not predict any (artificial) 
instability. The dependence of uj on A, obtained from Eq. (15), is shown by a solid line in Fig. 3. In the same figure, 
crosses show resonant values of the frequency found from numerical simulations of Eq. (3) with a small coefficient 
ai in front of the variable part of the nonlinearity coefficient, see Eq. (4). In the simulations, the forcing frequency 
u was varied at the fixed small a\, with the purpose to identify a value that generated strongest resonant response. 
The relative difference between the thus found resonance frequency and the value predicted by Eq. (15) is about 0.1, 
and the overall behavior of the curves is identical. It is worthy to note that ujq almost coincides with the soliton's 
frequency |%|. Thus, the results presented in Figs. 1 and 3 justify the validity of the VA based on the ansatz (9). 

The phase plane of Eqs. (12) and (13) for ao — 1 and a\ = is shown in Fig. 4(a), where arrows point out a direction 
of motion along a trajectory. The phase plane for the case of ao = — 1 is obtained by the shift b — > tt — 6, while that 
for the case ao = is obtained by setting a F p — > 0. As it follows from here, the stable FP, which corresponds to the 
discrete soliton, exists for either sign of a , and vanishes if ao = 0. As it also follows from Fig. 4(a), the evolution 
initiated by the initial condition with a(0) = app and small 6(0) corresponds to oscillations near the FP. However, 
for large values of |6(0)|, the asymptotic value of a(t) at t — > oo tends to zero. This fact is in qualitative agreement 
with the above-mentioned result that the addition of a chirp may destroy the soliton. 



C. The variational approximation for the nonstationary model 

We now proceed to the case of the ac-driven system, with a\ ^ 0. If a\ is small, strong response of the system 
to the time-periodic modulation is expected when the modulation frequency u> is close to the eigenfrequency lo$ of 
the internal oscillations of the soliton in the unperturbed system, which is given by Eq. (15); in fact, the resonant 
response was already taken into regard when collecting the data shown by crosses in Fig. 3. Moreover, the dynamics 
is expected to become chaotic, via the resonance-overlapping mechanism, if the modulation amplitude ai exceeds 
some threshold value. 
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The Poincare map illustrating a typical example of the chaotic behavior, as found from the numerical solution 
of Eqs. (12) and (13), is presented in Fig. 4(b). Shown are the discrete trajectories initiated by sets of the initial 
conditions, namely, the one with (b\,a\) = (0,0.789), that corresponds to the stationary discrete soliton with A = 
1 in the unperturbed system (ai = 0), and (62,^2) = (0.13,0.74). The modulation frequency u is close to the 
eigenfrequency of small oscillations uj - For the former initial condition, the point in the space (b, a) is chaotically 
moving away from the unperturbed FP. However, the chaotic evolution is a transient feature, as the point eventually 
moves so that a(t) asymptotically tends to zero, implying infinite broadening of the soliton. 

As for the second set of initial conditions, a new FP is found in a vicinity of the unperturbed one. This new FP 
predicts the existence of quasi-stationary discrete FRM solitons in the case of the slow modulation. Similar behavior 
near the corresponding stationary point is observed for the case ao = — 1. 

D. Direct simulations 

We have performed systematic comparison of the predictions produced by the VA against direct simulations of the 
full DNLS equation (3). The simulations show that, generally speaking, VA correctly predicts only an initial stage 
of the dynamics. The radiation of linear waves by a soliton, which is ignored by the VA, gives rise to an effective 
dissipation, that makes the resonance frequency different from luq. Furthermore, since luq depends on W, and the 
radiation loss results in gradual decrease of W, the soliton decouples from the resonance. In principle, VA might 
be made more accurate by adding a radiation mode ("tail") to the ansatz, cf. the analysis developed in Ref. 23 for 
the soliton in the continuum NLS equation (see also the review 14 ), but we do not aim to develop such an involved 
generalization of the VA in the present work. In any case, a conclusion is that the dynamics of the discrete soliton, as 
found from direct numerical simulation of Eq. (3) for a\ < 0.05, is close to that predicted by the variational equations 
(12) and (13). Namely, oscillations of the soliton's parameters are regular for very small modulations, and become 
chaotic when a\ exceed a threshold, see below. 

Typical examples of the soliton dynamics with u = 0.5 and different values of the modulation amplitude a\ are 
displayed in Fig. 5. An important observation, which is not predicted at all by the single-soliton ansatz, is splitting 
of the pulse, which is observed in Fig. 5(b). Note that for other values of a\, in Figs. 5(a) and 5(c), a stable soliton 
is observed, centered at n = 0, whose parameters oscillate because of the modulation. Therefore, the splitting which 
occurs at a\ > 0.1 is due to an interplay between the soliton itself, its intrinsic eigenmodes, and the energy exchange 
with radiation modes (continuous spectrum). It is noteworthy that the splitting is qualitatively similar to that 
revealed by direct simulations of the continuum NLS equation with a term accounting for periodic modulation of the 
linear dispersion term [whose discrete counterpart is the finite-difference combination in Eq. (3)], which was reported 
in Ref. 24 . A similar phenomenon was also observed in the discrete model with the finite-difference term subject to 
periodic modulation 3 . 

Results of the systematic numerical study of the splitting of the pulse with the initial amplitude A = 1 arc 
summarized in Fig. 6. In the simulations, absorbing boundary conditions, were used, the total number of particles 
was N > 200, and the dimensionless simulation time was, at least, 60tt/lo. We classify as splitting cases when at least 
two pulses emerge, moving in opposite directions, and no pulse with an appreciable amplitude stays around n = 0. 
For ai > 0.2, the modulation results in generation of several moving pulses. However, if a soliton with conspicuous 
amplitude is eventually found around n = 0, this case was classified as a "stable soliton" . 

Figure 6 also displays the dependence of a threshold amplitude a\ , past which the initial state chaotically drifts to 
a = 0, versus u> is also presented, as found from simulations of Eqs. (12) and (13). As is seen, the splitting actually 
occurs far above the threshold in a region of the developed dynamical chaos. The diagram for the case ao = — 1 looks 
similar, but not exactly the same. 

IV. THE AVERAGED EQUATION FOR THE CASE OF RAPID MODULATIONS 

In this section we consider the case of high-frequency modulations, with e = 1/uj <C 1. Note that we do not require 
ai to be small. In this case, it is natural to use the multiscale method 15 ' 11 . To this end, we introduce a set of time 
scales t = t/e, tk = £ k t, where k = 0, 1,2 . . ., and look for a solution in the form 

u n = U n + + e 2 u^ + ... (16) 
We substitute Eq. (16) into Eq. (3) and collect terms at the same order in e. Then, at order e° we obtain 
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a solution to which is 



+ + o (^"+1 + - 2(7 «) + 0Wl^n| 2 ^n - , (17) 

where a(r) = a(t/e), and [/„ is a function of the slow variables tk- After averaging in the fast variable r, one has 

*^ + ^ (C^n+i + C^n-i - 2C/„) + a |^| 2 C/„ = , (18) 

where a = (a(r)) standing for the average value of the variable coefficient a(r). Then the equation for first correction 
takes the form 

8 (1) 

*^T- = -[«(*) -«o] |(7 n | 2 t/„, 

l#> =i(/Xl-(Ml))|f/n| 2 f/n, 

where iti = /J"[a(x) — a^dx, and (...) again stands for the average value. 
At order e 1 , we obtain d\J n jdt\ = 0, and 

«£° = (/1 2 - (A*2» [|t/„| 2 ((7 n +l + 0„-l)- 

i[/2(^ +1 + L/^) - l|C/„ +1 | 2 C/„ +1 - ^C/^rl 2 ^-! 

-^[(Mi-(^i» 2 -2M]|C/ n | 4 ;7„ , 

where /i2 = J*J~ [/"l (#) — (ni)]dx, and M = ((a* 2 ) — (a*i) 2 )/2- Finally, at order e 2 we find 

= iM[\U n+1 \ 2 (2\U n \ 2 U n+1 + UlU* n+1 ) + 
\U n . 1 \ 2 {2\U n \ 2 U n . 1 + U 2 n U*_ 1 )- 

3\u n \\u n+1 + L/ n _i)] + 2 i Ma |;y„| 6 ;y„ (19) 

Substituting Eqs. (18) and (19) into the relation 

du n du n du n 2 du n 

e^r. 1- e 



9i dto dt\ dti 
one can derive the averaged equation, 

iUn + \{U n +i + Un-i - 2U n ) + a \U n \ 2 U n = 

-2Ma e 2 \U n \ 6 U n - Me 2 
x [\U n+1 \ 2 (2\U n \ 2 U n+1 + U 2 U* +1 ) + 
\U n - 1 \ 2 (2\U n \ 2 U n ^ 1 + U 2 U:_ 1 )- 

3\U n \ 4 U n+1 - 3\U n \ 4 U n ^] , (20) 

where M = a 2 / 4 for the case of the periodic modulation in Eq. (4). 

Equation (20) is the higher-order DNLS equation produced by the averaging procedure, which contains extra on- 
site and intersite (nonlocal) nonlinearities. A change of variables q n = U n + e 2 M\U n \ A U n allows to rewrite Eq. (20), 
retaining only terms up to 0(e 2 ), in the following form 

iQn + 7^(q n +i + q n -i - 2q n ) + a |<7„| 2 <7„ = 

ie 2 M [3\q n \\q n+1 + qn-!) + 2\q n \ 2 q 2 n {q* n+1 + q*^^ 

\q n +i\ 4 q n +i + \q n -i\ 4 qn-i] - 

e 2 M[\q n+1 \ 2 (2\q n \ 2 q n+1 + q n q* n+1 ) + 

\q n - 1 \ 2 (2\q n \ 2 q n ^+q 2 n q* n _ 1 )] . (21) 
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An advantage of the equation in the form (21) is that it can be derived from a Lagrangian, 



1 GO 

L q = L --e 2 M J2 (kn+i| 2 -k| 2 ) 2 



2 

n— — oo 

x{QnQn+i +q n qn+i) , (22) 

where L is obtained from the underlying Lagrangian (8) by the substitution u n — ► g„ and a(t) — > ao- The existence 
of the Lagrangian L g allows one to apply the variational approximation (VA) like in Section III. 
For the application of VA, we take the ansatz for q n in the form 

q n = Bexp(iip + ic\n\- /3\n\), (23) 

cf. Eq. (9). Substituting Eq. (23) into Eq. (22), we calculate the effective Lagrangian, 

sinh 2 (/3) tanh 3 (/3) 



L q = L - 4:€ MWq cos(c) 



sinh(3/3) 



Here Lq is the same expression as in Eq. (10), with a change b —> c, a — > f3, W — > W q = B 2 coth(f3), and a(t) — > ao. 
Now one can deduce a dynamical system for the variable c and /? similar to Eqs. (12) and (13). The fixed point 
(/?fpj 0) for a = 1, or (/3 FP , 7r) for a = —1, of this system represents a FRM soliton in the case of rapid modulations, 
where /3fp is to be found from the equation 

sinh(/? FP ) - ^[1 + 3tanh 2 (/? F p)] + 

4sign(a )e 2 MVF 2 sinh(/3 F p) tanh 2 (/3 F p) 

[10 + 15 cosh(2/3 F p) - cosh(4 /3 F p)] , , 

X [l + 2cosh(2/3 F p)] 2 1 j 

The norm W and amplitude A of the field U n in the averaged soliton are related to those of the field q n as 

Wt*W q [l- 2e 2 MW 2 (tanh 3 /?) coth(3/3)] , 



Ax> A q {l-e 2 MA A ). (25) 
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The dependence W(A) found from Eqs. (24) and (25) at different values of S = af/(4w 2 ) is displayed in Fig. 7. 
Different curves in the figure terminate at finite values of A because the relation (25), as well as the change of 
variables U n — ► q n , arc not valid outside the corresponding intervals. As it is suggested by Fig. 7, one can effectively 
control the soliton by an appropriate choice of the modulation parameters. Increase of the total number of particles 
in the averaged soliton, as compared to that in the unperturbed soliton with the same amplitude, is clearly seen in 
Fig. 7. 



V. CONCLUSIONS 



We have studied the dynamics of an array of Bosc-Einstein condensates with the time-dependent scattering length. 
Applying the variational approximation, the frequency lvq of small intrinsic oscillations of the soliton was predicted. 
The possibility of chaotic dynamics in the near-resonance case, when the driving frequency uj is close to uio, was shown. 
Direct simulations have demonstrated that the modulations of sufficient strength may result in splitting of the soliton. 
Results of the simulations were summarized in the form of the diagram which shows the splitting regions in the (w, a\) 
plane. The existence of stable Feshbach-resonance-managed discrete matter-wave solitons was demonstrated in the 
cases of both slow and rapid modulation of the nonlincarity coefficient. In the latter case, the soliton dynamics reduces 
to the generalized DNLS equation, which involves additional on-site and inter-site nonlinearities. By making use of 
this equation, properties of the averaged soliton were predicted. In particular, increase of the total number of atoms 
in this soliton in comparison with the ordinary discrete soliton of the same amplitude was shown. 

The chirp imprinting discussed in Sect. II can be an effective tool, similar to the phase-engineering method 25 , 
for manipulating the condensate's wave function. Pulse splitting induced by the chirp imprinting, or otherwise by 
the application of Fcshbach-resonance modulation, can be used as a source of coherent pulse pairs in an atomic 
Mach-Zchnder interferometer 26 . 
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FIG. 1. (a) The inverse width a (left axis) and the frequency \ (right axis) of the soliton vs. its amplitude A in the DNLS 
model without the time- modulation, do = 1. (b) The norm W (left axis) and the soliton's area S (right axis) vs. A. Point 
symbols represent data found from the numerical solution of the nonlinear eigenvalue problem; the solid lines are the prediction 
of the analytical variational approximation [see Eqs.(ll) and (14)]. 
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FIG. 2. The critical value of the chirp added to the fundamental discrete soliton, see Eq. (7), which splits the soliton into 
two separating pulses, vs. the amplitude of the unperturbed fundamental soliton. Squares (pluses) corresponds to oo = 1 
(oo = -1). 




FIG. 3. The frequency of small intrinsic oscillations of the discrete solitons around the stationary configurations, in the case 
ao = 1, vs. the soliton's amplitude A. The solid line shows the frequency too as predicted, in the framework of the variational 
approximation, by Eq. (15). Points connected by the dotted line are values of the forcing frequency which produce a resonant 
response in numerical simulations of Eq. (3) with a small time-periodic forcing term added to it. For comparison, the dashed 
line shows the soliton's internal frequency |x|- 
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FIG. 4. (a) The phase plane of the dynamical system based on Eqs. (12) and (13), in the case of ao = 1, oi = 0, and 
W = 1.5202. Such a value of W corresponds to a soliton with A = 1. (b) An example of chaotic dynamics for the periodi- 
cally-modulated system at W = 1.5202, a = 1, oi = 0.02766 and w = 0.481. 
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FIG. 5. Evolution of a discrete soliton with the initial amplitude A — 1 in the periodically modulated system with 
ao = l,w = 0.5, and different values of oi. 
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FIG. 6. The diagram in the plane (lo, ai), for the case ao = 1. Open and solid rectangles correspond to stable and splitting 
solitons, respectively. The initial soliton's amplitude is A — 1. The solid line is the chaos-onset threshold as predicted by the 
variational equations. 
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FIG. 7. The dependence of W vs. A of an average soliton (dashed lines) is compared with that of the unperturbed DNLS 
equation (solid line), ao = 1. 
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